#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 19 22:54:11 2023

@author: liqingsimac
"""

import numpy as np
import matplotlib.pyplot as plt

x=np.array([7.0, 10.5, 13.0, 17.5, 34.0, 40.5, 44.5, 48.0, 56.0,
            61.0, 68.5, 76.5, 80.5, 91.0, 96.0, 101.0, 104.0, 106.5, 
            111.5, 118.0, 123.5, 136.5, 142.0, 146.0, 150.0, 157.0, 158.0])

y1=np.array([44, 45, 47, 50, 50, 38, 30, 30, 34, 
             36, 34, 41, 45, 46, 43, 37, 33, 28, 
             32, 65, 55, 54, 52, 50, 66, 66, 68])

y2=np.array([44, 59, 70, 72, 93, 100, 110, 110, 110,
             117, 118, 116, 118, 118, 121, 124, 121, 121, 
             121, 122, 116, 83, 81, 82, 86, 85, 68])

##画出点、连线
plt.plot(x,y1,'.-',label='y=f(x)')
plt.plot(x,y2,'.-',label='y=g(x)')
plt.legend(loc='upper left')

##画出点、连线、竖线
plt.plot(x,y1,'.-')
plt.plot(x,y2,'.-')
plt.plot([x,x],[y1,y2])

##使用循环，计算一个个的梯形面积
N=len(x)
Area=0
for k in range(N-1):
    LeftSide = y2[k]-y1[k]
    RightSide = y2[k+1]-y1[k+1]
    dx = x[k+1]-x[k]
    Area=Area + (LeftSide + RightSide)*dx/2

##使用数组计算
VB=y2-y1
Area=np.sum((VB[:-1]+VB[1:])*np.diff(x)/2)

RealArea=Area*(40/18)**2
print('The area by linear interpolation is %.2f square kilometers.'%RealArea)

##三次样条插值
from scipy.interpolate import interp1d
a=np.min(x)
b=np.max(x)
n=len(x)*5
xnew=np.linspace(a,b,n)
f1new=interp1d(x,y1,'cubic')
y1new=f1new(xnew)
f2new=interp1d(x,y2,'cubic')
y2new=f2new(xnew)

##画出三次样条插值得到的两个函数图像
plt.plot(x,y1,'.')
plt.plot(x,y2,'.')
plt.plot(xnew,y1new,'-')
plt.plot(xnew,y2new,'-')

##使用三次样条插值数据计算面积
VertBars=y2new-y1new
SpArea=np.sum((VertBars[:-1]+VertBars[1:])*np.diff(xnew)/2)

RealSpArea=SpArea*(40/18)**2
print('The area by spline interpolation is %.2f square kilometers.'%RealSpArea)

